02-2022

Bonjour!

Aujourd’hui

  • Analyses statistiques de base avec R
  • Un peu de révision au passage: dplyr, ggplot2
  • Très peu de théorie! On assume que vous avez fait un cours universitaire de 1er cycle en statistiques et que vous êtes familier.e avec les termes “test de student”, “ANOVA”, “régression”…
  • Masteriser la fonction lm()
  • Code disponible ici:
  • Défis!

Données pour l’atelier

Nous allons réutiliser les ensembles de données que vous avez utilisé lors de deux atelier précédents: gapminder, palmer pinguins, starwars. Nous devons également charger les librairies dplyr et ggplot2.

library(dplyr)
library(ggplot2)
library(palmerpenguins)
library(gapminder)
data("starwars")
data("penguins")
data("gapminder")
# #au besoin
# install.packages(c("penguins","gapminder"))

Types de variables

  • Selon StatCan
  • variables catégoriques/qualitatives: nominales & ordinales
  • variables numériques/quantitatives: continues & discrètes
head(penguins)
## # A tibble: 6 × 8
##   species island bill_length_mm bill_depth_mm flipper_length_… body_mass_g sex  
##   <fct>   <fct>           <dbl>         <dbl>            <int>       <int> <fct>
## 1 Adelie  Torge…           39.1          18.7              181        3750 male 
## 2 Adelie  Torge…           39.5          17.4              186        3800 fema…
## 3 Adelie  Torge…           40.3          18                195        3250 fema…
## 4 Adelie  Torge…           NA            NA                 NA          NA <NA> 
## 5 Adelie  Torge…           36.7          19.3              193        3450 fema…
## 6 Adelie  Torge…           39.3          20.6              190        3650 male 
## # … with 1 more variable: year <int>

Défi!

Identifiez les variables qualitatives et quantitatives de l’ensemble de données starwars

Défi!

Identifiez les variables qualitatives et quantitatives de l’ensemble de données starwars

head(starwars)
glimpse(starwars)
## Rows: 87
## Columns: 14
## $ name       <chr> "Luke Skywalker", "C-3PO", "R2-D2", "Darth Vader", "Leia Or…
## $ height     <int> 172, 167, 96, 202, 150, 178, 165, 97, 183, 182, 188, 180, 2…
## $ mass       <dbl> 77.0, 75.0, 32.0, 136.0, 49.0, 120.0, 75.0, 32.0, 84.0, 77.…
## $ hair_color <chr> "blond", NA, NA, "none", "brown", "brown, grey", "brown", N…
## $ skin_color <chr> "fair", "gold", "white, blue", "white", "light", "light", "…
## $ eye_color  <chr> "blue", "yellow", "red", "yellow", "brown", "blue", "blue",…
## $ birth_year <dbl> 19.0, 112.0, 33.0, 41.9, 19.0, 52.0, 47.0, NA, 24.0, 57.0, …
## $ sex        <chr> "male", "none", "none", "male", "female", "male", "female",…
## $ gender     <chr> "masculine", "masculine", "masculine", "masculine", "femini…
## $ homeworld  <chr> "Tatooine", "Tatooine", "Naboo", "Tatooine", "Alderaan", "T…
## $ species    <chr> "Human", "Droid", "Droid", "Human", "Human", "Human", "Huma…
## $ films      <list> <"The Empire Strikes Back", "Revenge of the Sith", "Return…
## $ vehicles   <list> <"Snowspeeder", "Imperial Speeder Bike">, <>, <>, <>, "Imp…
## $ starships  <list> <"X-wing", "Imperial shuttle">, <>, <>, "TIE Advanced x1",…

Le modèle linéaire général

  • Une seule variable dépendante (réponse) quantitative
  • Distribution normale (plus ou moins)

\[ y = mx + b \] \[ y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_px_{ip} + \epsilon_i \] \[ \epsilon \sim N(0,\sigma^2) \]

Le modèle linéaire général

  • Une seule variable dépendante (réponse) quantitative
  • Distribution normale (plus ou moins)
hist(rnorm(10000))

Le modèle linéaire général

Test de Student (test de t)

  • Différence entre deux moyennes
  • Chez le manchot Adélie, est-ce que la masse corporelle diffère entre les mâles et les femelles?
penguins %>%
  filter(species == 'Adelie') %>% 
  group_by(sex) %>% 
  summarize(masse.moy = mean(body_mass_g))
## # A tibble: 3 × 2
##   sex    masse.moy
##   <fct>      <dbl>
## 1 female     3369.
## 2 male       4043.
## 3 <NA>         NA

Test de Student (test de t)

penguins %>%
  filter(species == 'Adelie', sex != 'NA') %>%
  ggplot(aes(x=body_mass_g, colour=sex)) + geom_density()

Test de Student (test de t)

Il semble y avoir une différence de masse entre les sexes. La fonction t.test nous permet de valider cette hypothèse.

adelie <- penguins %>%
  filter(species == 'Adelie', sex != 'NA')
t.test(body_mass_g ~ sex, data = adelie)
## 
##  Welch Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -13.126, df = 135.69, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -776.3012 -573.0139
## sample estimates:
## mean in group female   mean in group male 
##             3368.836             4043.493

Test de Student (test de t)

Il semble y avoir une différence de masse entre les sexes. La fonction t.test nous permet de valider cette hypothèse.

adelie <- penguins %>%
  filter(species == 'Adelie', sex != 'NA')
t.test(body_mass_g ~ sex, data = adelie)
## 
##  Welch Two Sample t-test
## 
## data:  body_mass_g by sex
## t = -13.126, df = 135.69, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -776.3012 -573.0139
## sample estimates:
## mean in group female   mean in group male 
##             3368.836             4043.493
format(2.2e-16, scientific = F)
## [1] "0.00000000000000022"

Test de Student (test de t)

La fonction lm fonctionne pour tous les types de modèles linéaires (incluant le test de t).

modele <- lm(body_mass_g ~ sex, data = adelie)
summary(modele)
## 
## Call:
## lm(formula = body_mass_g ~ sex, data = adelie)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -718.49 -218.84  -18.84  225.00  731.51 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3368.84      36.34   92.69   <2e-16 ***
## sexmale       674.66      51.40   13.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 310.5 on 144 degrees of freedom
## Multiple R-squared:  0.5447, Adjusted R-squared:  0.5416 
## F-statistic: 172.3 on 1 and 144 DF,  p-value: < 2.2e-16

Test de Student (test de t)

\[ Y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \] Quelles valeurs la variable “x” (sex) peut-elle prendre?

summary(modele)
## 
## Call:
## lm(formula = body_mass_g ~ sex, data = adelie)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -718.49 -218.84  -18.84  225.00  731.51 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3368.84      36.34   92.69   <2e-16 ***
## sexmale       674.66      51.40   13.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 310.5 on 144 degrees of freedom
## Multiple R-squared:  0.5447, Adjusted R-squared:  0.5416 
## F-statistic: 172.3 on 1 and 144 DF,  p-value: < 2.2e-16

Défi!

Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo? On peut extraire les données pertinentes comme suit:

nab_tat <- starwars %>% filter(homeworld %in% c('Naboo','Tatooine'))
ggplot(aes(x=homeworld,y=height),data=nab_tat) + geom_boxplot()

Défi!

Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo?

t.test(height ~ homeworld, data = nab_tat)
## 
##  Welch Two Sample t-test
## 
## data:  height by homeworld
## t = 0.42386, df = 18.947, p-value = 0.6764
## alternative hypothesis: true difference in means between group Naboo and group Tatooine is not equal to 0
## 95 percent confidence interval:
##  -22.27276  33.58185
## sample estimates:
##    mean in group Naboo mean in group Tatooine 
##               175.4545               169.8000

Défi!

Dans les données starwars, est-ce que les personnages provenant de Tatooine sont plus grands que ceux provenant de Naboo?

modele <- lm(height ~ homeworld, data = nab_tat)
summary(modele)
## 
## Call:
## lm(formula = height ~ homeworld, data = nab_tat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.455  -6.800   7.545  13.200  48.545 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        175.455      9.276  18.914 8.76e-14 ***
## homeworldTatooine   -5.655     13.443  -0.421    0.679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.77 on 19 degrees of freedom
## Multiple R-squared:  0.009227,   Adjusted R-squared:  -0.04292 
## F-statistic: 0.1769 on 1 and 19 DF,  p-value: 0.6787

ANOVA

Lorsque notre variable qualitative indépendante comprend plus de deux niveaux, nous pouvons utiliser une analyse de variance (ANOVA) plutôt qu’un test de t. On ajuste le modèle de la même façon avec la fonction lm:

modele <- lm(body_mass_g ~ species, data = penguins)
summary(modele)
## 
## Call:
## lm(formula = body_mass_g ~ species, data = penguins)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1126.02  -333.09   -33.09   316.91  1223.98 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3700.66      37.62   98.37   <2e-16 ***
## speciesChinstrap    32.43      67.51    0.48    0.631    
## speciesGentoo     1375.35      56.15   24.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 462.3 on 339 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.6697, Adjusted R-squared:  0.6677 
## F-statistic: 343.6 on 2 and 339 DF,  p-value: < 2.2e-16

ANOVA

Les coefficients (béta) du modèle indiquent la différence de masse entre chaque niveau et le niveau de référence. Par défaut, le niveau de référence d’un facteur est celui qui arrive en première position lorsqu’on classe les niveaux en ordre alphabétique. Pour nos manchots, il y a trois espèces: Adelie, Chinstrap, Gentoo. Adelie est donc le niveau de référence.

levels(penguins$species)
## [1] "Adelie"    "Chinstrap" "Gentoo"
coef(modele)
##      (Intercept) speciesChinstrap    speciesGentoo 
##       3700.66225         32.42598       1375.35401

masse moyenne Adelie = 3700, Chinstrap = 3700+32, Gentoo = 3700+1375

ANOVA

masse moyenne Adelie = 3700, Chinstrap = 3700+32, Gentoo = 3700+1375

ggplot(aes(x=species,y=body_mass_g),data=penguins) + geom_boxplot()

ANOVA

Souvent, pour une ANOVA, on cherche plutôt à tester si la variable indépendante influence au moins un des niveaux, peu importe lequel, et sans nécessairement comparer les niveaux à un niveau de référence. On veut un tableau d’ANOVA qui ressemble à cela:

https://online.stat.psu.edu/stat415/lesson/13/13.2

ANOVA

Il suffit d’utiliser la fonction anova pour obtenir le tableau d’ANOVA de notre modèle:

anova(modele)
## Analysis of Variance Table
## 
## Response: body_mass_g
##            Df    Sum Sq  Mean Sq F value    Pr(>F)    
## species     2 146864214 73432107  343.63 < 2.2e-16 ***
## Residuals 339  72443483   213698                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA

On complémente souvent le test de F de l’ANOVA avec des comparaisons par paires, pour déterminer quels groupes diffèrent les uns des autres. Dans ce cas, y a-t-il une différence entre Chinstrap et Adelie?

ggplot(aes(x=species,y=body_mass_g),data=penguins) + geom_boxplot()

ANOVA

Comparaisons par paires à l’aide de tests de Tukey (“Tukey’s Honest Significant Difference”):

TukeyHSD(aov(modele))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = modele)
## 
## $species
##                        diff       lwr       upr     p adj
## Chinstrap-Adelie   32.42598 -126.5002  191.3522 0.8806666
## Gentoo-Adelie    1375.35401 1243.1786 1507.5294 0.0000000
## Gentoo-Chinstrap 1342.92802 1178.4810 1507.3750 0.0000000

Défi!

À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?

donnees_2007 <- gapminder %>% filter(year == 2007)
ggplot(aes(x=continent,y=lifeExp),data=donnees_2007) + geom_boxplot()

Défi!

À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?

modele <- lm(lifeExp ~ continent, donnees_2007)
anova(modele)
## Analysis of Variance Table
## 
## Response: lifeExp
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## continent   4 13060.7  3265.2  59.714 < 2.2e-16 ***
## Residuals 137  7491.2    54.7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Défi!

À votre tour! Utilisez les données gapminder pour ajuster un modèle d’ANOVA. En 2007, y avait-il une différence d’espérance de vie (lifeExp) entre les continents (continent)? Le cas échéant, quelle(s) paire(s) de continents montraient une différence pour cette variable?

TukeyHSD(aov(modele))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = modele)
## 
## $continent
##                       diff        lwr       upr     p adj
## Americas-Africa  18.802082  13.827042 23.777121 0.0000000
## Asia-Africa      15.922446  11.372842 20.472051 0.0000000
## Europe-Africa    22.842562  18.155858 27.529265 0.0000000
## Oceania-Africa   25.913462  11.183451 40.643472 0.0000305
## Asia-Americas    -2.879635  -8.299767  2.540497 0.5844901
## Europe-Americas   4.040480  -1.495233  9.576193 0.2630707
## Oceania-Americas  7.111380  -7.910342 22.133102 0.6862848
## Europe-Asia       6.920115   1.763372 12.076858 0.0027416
## Oceania-Asia      9.991015  -4.895221 24.877251 0.3464970
## Oceania-Europe    3.070900 -11.857807 17.999607 0.9793776

Défi!

Merveilleux, mais… en regardant ces boîtes à moustaches (boxplot), trouvez-vous qu’il y a quelque chose de louche?

ggplot(aes(x=continent,y=lifeExp),data=donnees_2007) + geom_boxplot()

Suppositions des modèles linéaires

\[ Y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \]

\[ \epsilon \sim N(0,\sigma^2) \]

Quelques suppositions des modèles linéaires:

  • Les points de données sont indépendants (p. ex., on ne mesure pas le même manchot ou le même pays 10 fois!)
  • Les résidus du modèle suivent une distribution normale
  • La variance des résidus est identique pour tous les groupes (homoscédasticité)

Suppositions des modèles linéaires

Nous pouvons vérifier ces suppositions à l’aide de tests. IMPORTANT: pour ces deux tests, nous voulons que p > 0.05!

#1. normalité des résidus
mod.resid <- resid(modele)  # Sortir les résidus
shapiro.test(mod.resid)  # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## 
##  Shapiro-Wilk normality test
## 
## data:  mod.resid
## W = 0.96708, p-value = 0.001696

Suppositions des modèles linéaires

Nous pouvons vérifier ces suppositions à l’aide de tests. IMPORTANT: pour ces deux tests, nous voulons que p > 0.05!

#1. normalité des résidus
mod.resid <- resid(modele)  # Sortir les résidus
shapiro.test(mod.resid)  # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## 
##  Shapiro-Wilk normality test
## 
## data:  mod.resid
## W = 0.96708, p-value = 0.001696
# 2. homoscédasticité
bartlett.test(lifeExp ~ continent, donnees_2007) # Test de Bartlett: p < 0.05 = MAUVAIS!
## 
##  Bartlett test of homogeneity of variances
## 
## data:  lifeExp by continent
## Bartlett's K-squared = 45.85, df = 4, p-value = 2.646e-09

Suppositions des modèles linéaires

Il faut vérifier les suppositions d’un modèle avant de l’interpréter! S’il y a un problème, on ne regarde même pas les résultats.

Solutions potentielles:

  • exclure le(s) groupe(s) particulièrement problématique(s), ce qui diminue la portée de l’analyse
  • transformer la variable réponse (p.ex. log)
  • changer la structure du modèle en ajoutant des variables indépendantes
  • utiliser un modèle linéaire généralisé ou d’autres outils statistiques plus avancés

ANOVA avec deux facteurs

Une ANOVA peut inclure plus d’une variables qualitatives, ce qui permet notamment de quantifier les intéractions entre deux ou plusieurs variables. Par exemple, est-ce que l’effet du sexe sur la masse des manchots dépend de l’espèce?

penguins <- tidyr::drop_na(penguins)
ggplot(aes(x=species,y=body_mass_g,fill=sex),data=penguins) + geom_boxplot()

ANOVA avec deux facteurs

Une ANOVA peut inclure plus d’une variables qualitatives, ce qui permet notamment de quantifier les intéractions entre deux ou plusieurs variables. Par exemple, est-ce que l’effet du sexe sur la masse des manchots dépend de l’espèce?

modele <- lm(body_mass_g ~ species + sex, data = penguins)
anova(modele)
## Analysis of Variance Table
## 
## Response: body_mass_g
##            Df    Sum Sq  Mean Sq F value    Pr(>F)    
## species     2 145190219 72595110  724.21 < 2.2e-16 ***
## sex         1  37090262 37090262  370.01 < 2.2e-16 ***
## Residuals 329  32979185   100241                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Que manque-t-il?

ANOVA avec deux facteurs

modele <- lm(body_mass_g ~ species * sex, data = penguins)
anova(modele)
## Analysis of Variance Table
## 
## Response: body_mass_g
##              Df    Sum Sq  Mean Sq F value    Pr(>F)    
## species       2 145190219 72595110 758.358 < 2.2e-16 ***
## sex           1  37090262 37090262 387.460 < 2.2e-16 ***
## species:sex   2   1676557   838278   8.757 0.0001973 ***
## Residuals   327  31302628    95727                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA avec deux facteurs

N’oublions pas qu’il s’agit encore du modèle linéaire général:

\[ Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_px_{ip} + \epsilon_i \]

summary(modele)
## 
## Call:
## lm(formula = body_mass_g ~ species * sex, data = penguins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -827.21 -213.97   11.03  206.51  861.03 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               3368.84      36.21  93.030  < 2e-16 ***
## speciesChinstrap           158.37      64.24   2.465  0.01420 *  
## speciesGentoo             1310.91      54.42  24.088  < 2e-16 ***
## sexmale                    674.66      51.21  13.174  < 2e-16 ***
## speciesChinstrap:sexmale  -262.89      90.85  -2.894  0.00406 ** 
## speciesGentoo:sexmale      130.44      76.44   1.706  0.08886 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 309.4 on 327 degrees of freedom
## Multiple R-squared:  0.8546, Adjusted R-squared:  0.8524 
## F-statistic: 384.3 on 5 and 327 DF,  p-value: < 2.2e-16

Régression & ANCOVA

Régression simple

Chez les personnages de star wars, y a-t-il une relation entre la masse corporelle et la grandeur?

ggplot(aes(x=height,y=mass),data=starwars) + geom_point()

Régression simple

Chez les personnages de star wars, y a-t-il une relation entre la masse corporelle et la grandeur?

ggplot(aes(x=height,y=mass),data=starwars) + geom_point()

Régression simple

Chez les personnages de star wars (autres que Jabba the Hut), y a-t-il une relation entre la masse corporelle et la grandeur?

sw <- starwars %>% filter(mass < 1000)
ggplot(aes(x=height,y=mass),data=sw) + geom_point()

Régression simple

Chez les personnages de star wars (autres que Jabba the Hut), y a-t-il une relation entre la masse corporelle et la grandeur?

modele <- lm(mass ~ height, data = sw)
summary(modele)
## 
## Call:
## lm(formula = mass ~ height, data = sw)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.382  -8.212   0.211   3.846  57.327 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -32.54076   12.56053  -2.591   0.0122 *  
## height        0.62136    0.07073   8.785 4.02e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.14 on 56 degrees of freedom
## Multiple R-squared:  0.5795, Adjusted R-squared:  0.572 
## F-statistic: 77.18 on 1 and 56 DF,  p-value: 4.018e-12

Régression simple

Validons les suppositions. Premièrement, la normalité des résidus:

mod.resid <- resid(modele)
shapiro.test(mod.resid)  # Test de Shapiro-Wilk: p < 0.05 = MAUVAIS!
## 
##  Shapiro-Wilk normality test
## 
## data:  mod.resid
## W = 0.912, p-value = 0.0004688

Ça ne passe pas…

Régression simple

On ré-essaie en transformant notre variable dépendante.

sw$log_masse <- log10(sw$mass)
modele <- lm(log_masse ~ height, data = sw)
mod.resid <- resid(modele)
shapiro.test(mod.resid) #yé!
## 
##  Shapiro-Wilk normality test
## 
## data:  mod.resid
## W = 0.96702, p-value = 0.1156

Régression simple

Pour l’homogénéité de la variance, il suffit de regarder la distribution des résidus en fonction des valeurs de y estimées par le modèle.

Graphique tiré de ce site web

Régression simple

Pour l’homogénéité de la variance, il suffit de regarder la distribution des résidus en fonction des valeurs de y estimées par le modèle. Ici, on espère observer un nuage de point sans aucune forme triangulaire ou tendance directionnelle évidente.

plot(modele, which = 1)

Régression simple

Tout est bien. On peut interpréter le modèle…

summary(modele)
## 
## Call:
## lm(formula = log_masse ~ height, data = sw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25043 -0.06977  0.01734  0.05502  0.22060 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.9911872  0.0690142   14.36   <2e-16 ***
## height      0.0048730  0.0003886   12.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1052 on 56 degrees of freedom
## Multiple R-squared:  0.7374, Adjusted R-squared:  0.7327 
## F-statistic: 157.2 on 1 and 56 DF,  p-value: < 2.2e-16

Régression simple

Tout est bien. On peut interpréter le modèle et le visualiser.

ggplot(aes(x=height,y=log_masse),data=sw) +
  geom_point() + 
  geom_smooth(method='lm')

Régression simple

\[ y_i = \beta_0 + \beta_1x_{i} + \epsilon_i \]

Plutôt que d’utiliser ggplot pour ajuster le modèle à nouveau, nous pouvons aussi extraire les coéfficients (béta) de notre modèle et les utiliser pour “dessiner” la droite.

coef(modele)
## (Intercept)      height 
## 0.991187162 0.004872992
b0 <-  coef(modele)[1]
b1 <-  coef(modele)[2]
c(b0,b1)
## (Intercept)      height 
## 0.991187162 0.004872992

Quelle est l’équation de notre modèle linéaire?

Régression simple

\[ y_i = 0.9912 + 0.0049 * height_i + \epsilon_i \]

ggplot(aes(x=height,y=log_masse),data=sw) +
  geom_point() + 
  geom_abline(intercept = b0, slope = b1)

Défi!

À votre tour! Dans les données penguins, y a-t-il une relation entre la longueur du bec (bill_length_mm) et la longueur des nageoires (flipper_length_mm)?

Défi!

modele <- lm(bill_length_mm ~ flipper_length_mm, penguins)
par(mfrow=c(1,2))
plot(modele, which = 1:2)

Défi!

summary(modele)
## 
## Call:
## lm(formula = bill_length_mm ~ flipper_length_mm, data = penguins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.6367 -2.6981 -0.5788  2.0663 19.0953 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -7.21856    3.27175  -2.206    0.028 *  
## flipper_length_mm  0.25482    0.01624  15.691   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.148 on 331 degrees of freedom
## Multiple R-squared:  0.4265, Adjusted R-squared:  0.4248 
## F-statistic: 246.2 on 1 and 331 DF,  p-value: < 2.2e-16

longueur du bec = -7.21856 + 0.25482 * longueur de la nageoire

Défi!

ggplot(aes(x=flipper_length_mm,y=bill_length_mm),data=penguins) +
  geom_point() + 
  geom_smooth(method='lm')

ANCOVA

Analyse de covariance: 1 variable qualitative, 1 variable quantitative.

ggplot(aes( x=body_mass_g, y=bill_length_mm, colour = species), data = penguins) +
  geom_point() + geom_smooth(method = "lm", aes(fill = species))

ANCOVA

ANCOVA

modele.sans.interaction <- lm(bill_length_mm ~ body_mass_g + species, data = penguins)
summary(modele.sans.interaction) #equation pour le manchot Gentoo?
## 
## Call:
## lm(formula = bill_length_mm ~ body_mass_g + species, data = penguins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8291 -1.6728  0.1244  1.5318  9.2904 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      24.908763   1.089730  22.858  < 2e-16 ***
## body_mass_g       0.003755   0.000289  12.991  < 2e-16 ***
## speciesChinstrap  9.908762   0.355289  27.889  < 2e-16 ***
## speciesGentoo     3.539179   0.499814   7.081 8.71e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.419 on 329 degrees of freedom
## Multiple R-squared:  0.806,  Adjusted R-squared:  0.8043 
## F-statistic: 455.8 on 3 and 329 DF,  p-value: < 2.2e-16

ANCOVA

modele.avec.interaction <- lm(bill_length_mm ~ body_mass_g * species, data = penguins)
summary(modele.avec.interaction) #equation pour le manchot Gentoo?
## 
## Call:
## lm(formula = bill_length_mm ~ body_mass_g * species, data = penguins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4164 -1.6630  0.0683  1.5444  9.3138 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  27.1129056  1.6324193  16.609  < 2e-16 ***
## body_mass_g                   0.0031599  0.0004371   7.228 3.48e-12 ***
## speciesChinstrap              5.0612873  3.3101846   1.529    0.127    
## speciesGentoo                -0.5750259  2.7941191  -0.206    0.837    
## body_mass_g:speciesChinstrap  0.0013028  0.0008832   1.475    0.141    
## body_mass_g:speciesGentoo     0.0009698  0.0006225   1.558    0.120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.414 on 327 degrees of freedom
## Multiple R-squared:  0.8081, Adjusted R-squared:  0.8051 
## F-statistic: 275.3 on 5 and 327 DF,  p-value: < 2.2e-16

Régression multiple

Régression multiple

En ajoutant une troisième variable à notre modèle, on arrive à un modèle de régression multiple. Attention: les variables indépendantes ne doivent pas être trop corrélées les unes aux autres (problème de la multicolinéarité).

modele <- lm(bill_length_mm ~ body_mass_g + species + year, data = penguins)
summary(modele)
## 
## Call:
## lm(formula = bill_length_mm ~ body_mass_g + species + year, data = penguins)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4997 -1.6727  0.0599  1.4895  9.5913 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -5.978e+02  3.270e+02  -1.828   0.0685 .  
## body_mass_g       3.751e-03  2.879e-04  13.031  < 2e-16 ***
## speciesChinstrap  9.935e+00  3.541e-01  28.053  < 2e-16 ***
## speciesGentoo     3.540e+00  4.978e-01   7.110 7.28e-12 ***
## year              3.101e-01  1.629e-01   1.904   0.0578 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.41 on 328 degrees of freedom
## Multiple R-squared:  0.8082, Adjusted R-squared:  0.8058 
## F-statistic: 345.5 on 4 and 328 DF,  p-value: < 2.2e-16

Régression multiple

On valide le modèle de la même façon.

par(mfrow=c(1,2))
plot(modele, which = 1:2)

Régression multiple

La représentation graphique d’un modèle de régression multiple est plus difficile. Souvent, nous nous contentons de montrer les coefficients du modèle sur un forest plot.

sjPlot::plot_model(modele)

Régression multiple

La représentation graphique d’un modèle de régression multiple est plus difficile. Souvent, nous nous contentons de montrer les coefficients du modèle sur un forest plot. Pour pouvoir comparer la taille des effets, il faut centrer-réduire toutes les variables quantitatives:

modele <- lm(bill_length_mm ~ scale(body_mass_g) + species + scale(year), data = penguins)
sjPlot::plot_model(modele)

Au delà du modéle linéaire général

  • La variable réponse ne suit pas une distribution normale: modèles linéaire généralisés
  • Les points de données ne sont pas indépendants: modèles mixtes
  • La relation entre x et y n’est pas linéaire: modèles additifs (GAM)
  • Plusieurs variables dépendantes: analyse multidimensionnelle

Quelques ressources utiles: